Lattice Boltzmann simulations of apparent slip and contact angle in hydrophobic 

micro-channels 
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In this paper, we applied the Shan-Chen multiphase Lattice Boltzmann method to simulate two 
different parameters, contact angle (a static parameter) and slip length (a dynamic parameter), 
and we proposed a relationship between them by fitting those numerical simulation results. By 
changing the values of the strength of interaction between fluid particles (G) and the strength 
of interaction between fluid and solid surface (Gads), we simulated a series of contact angles and 
slip lengths. Our numerical simulation results show that both G and Gads have little effects on 
the relationship between contact angle and slip length. Using the proposed relationship between 
slip length and contact angle, we further derived an equation to determine the upper limit of nano- 
particles' diameter under which drag-reduction can be achieved when using nano-particles adsorbing 
method. 



I. INTRODUCTION 

Contact angle is a static parameter of measuring the 
wettability of a liquid on a solid surface, and it can be 
easily measured. Slip length is a dynamic parameter of 
quantifying the non-zero velocity boundary condition of a 
liquid flowing over a solid surface. Determination of slip 
length is very important for calculation of drag and other 
hydrodynamic behaviors of fluid flowing through micro- 
channels or over nano-scale patterned surfaces. However, 
it is very difficult to directly measure the apparent slip 
length. 

Non-slip boundary condition is extremely successful in 
describing macro-scale viscous flows in engineering ap- 
plications for more than one hundred years [l|. This as- 
sumption is supported only by macroscopic experimen- 
tal results. However, a series of experiments [2|-|j] and 
numerical simulation results 0-0| indicate that this as- 
sumption does not hold at micro- and nano-scale, and 
a slip boundary condition should be applied. Based 
on this slip boundary effect, artificial super-hydrophobic 
surfaces have been widely used in industrial production 
and people's daily lives, for example, self-cleaning paints, 
roof tiles, fabrics and glass windows that can be cleaned 
by a simple rainfallQ and the nano-particles adsorbing 
method [2| in improving oil recovery are all in practice. 

In 1823, Navier proposed a slip boundary condition 
that the fluid velocity at a point on a surface is pro- 
portional to the shear rate at the same point, v(xt) — 
Sdv(x)/dx, where S is defined as the slip length, but its 
value is hard to determine. Molecular dynamics simu- 
lations (MDS) have been widely used to study the rela- 
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tionship between fluid slip and the properties of fluid and 
solid, and the resultsproved that there is boundary slip 
on microscopic scaleQ. However, the conventional MDS 
method has difficulty to determine the small flow veloc- 
ity because the nonlinear coupling of the small bulk flow 
velocity with the large peculiar velocity of the thermal 
motion 10] and has difficulty to simulate large size sys- 
tems [ll|. Chen et al investigated the Couette geometry 
flows by means of a two-phase mesoscopic Lattice Boltz- 
mann(LB) model, and the results show that there is a 
strong relationship between the magnitude of slip and 
the solid-fluid inter action [l2j], but it is a very difficult or 
even an impossible task to compute the exact slip in de- 
pendence of interaction and it is also difficult to apply to 
engineering. 

In this present work, we focus on investigating the 
effects of wall wettabilities on the slip length. With a 
general bounce-back no-slip boundary condition applied 
to the interface between fluids and solid surfaces, to- 
gether with the Shan-Chen multiphase model[l3[, the LB 
method is used to simulate the Poiseuille flow. The sim- 
ulation results reveal that the wetting properties of the 
wall have an important effect on the flow. We proposed 
a relationship between slip length and contact angle ac- 
cording to these numerical simulation results. Using this 
relationship, it is easy to estimate the slip length because 
the contact angle is a parameter that can be easily mea- 
sured. 

Nano-partical adsorption has been proved an effective 
drag-reduction method that can be widely used for en- 
hancing oil recovery and many other practical applica- 
tions. However, it is still not clear how to properly select 
the size or diameter of nano-particles. The final section 
considers drag reduction by hydrophobic nanoparticles 
(HNPs). Placing these particles on the walls of a chan- 
nel constricts the channel but increases the slip length. 
An overall reduction in flow resistance occurs if the latter 
effect outweighs the former. Using the relationship be- 
tween slip length and contact angle, we further derived an 
equation to determine the upper limit of nano-particles' 
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diameter under which drag-reduction can be achieved. 



II. NUMERICAL MODEL 

A. The LBGK model 

The LB method, which involves a single relaxation 
time in the Bhatnagar-Gross-Krook (BGK) collision 
operator 14j , is used here. The time evolution of this 
model can be written as 



/ t (x + Ci Ai, t + At) - / 4 (x, t) = -i[/,(x, t) - /Hx, *)], 

T 

(1) 

where /i(x, t) is the single-particle distribution function 
for fluid particles that move in the direction at (x, t), 
f^ q (x, t) is the equilibrium distribution function, At is 
time step of simulation, and parameter r is the relax- 
ation time characterizing the collision processes by which 
the distribution functions relax towards their equilibrium 
distributions. In the two-Dimensional (2D) nine-particle 
model, the equilibrium distribution function, (x, t), 
depends only on local density and velocity and can be 
chosen as the following form 
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c = Ax/ At is lattice velocity, Ax is lattice distance, and 
At is time step of simulation. 

The mass density p and momentum density pu of the 
fluid are calculated from the first and second moments of 
the distribution function, i.e., 
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And the relaxation time tunes the kinematic viscosity as 

1, 



The non-slip boundary condition at solid-fluid interfaces 
is realized by the bounce-back rule, where the momen- 
tum of the fluid particles that meet a solid wall is simply 
reversed. The bounce-back rule is simple and compu- 
tationally efficient, and enables fluid flow simulations in 
complex geometries. On the inlet and outlet boundaries, 
the periodic boundary conditions are used. 



B. Shan-Chen multiphase model 

To simulate non-ideal multiphase fluids, the attrac- 
tive or repulsive interaction among molecules, which is 
referred to as the non-ideal interaction, should be in- 
cluded in the LB model. There are many approaches 
to incorporate non-ideal interactions, such as color-fluid 
model, interparticle-potential model, free-energy model, 
mean-field theory model and so on. The interparticle- 
potential model proposed by Shan & Chen[TJ|, [l5[ is to 
mimic microscopic interaction forces between the fluid 
components. This model modified the collision operator 
by using an equilibrium velocity that includes an inter- 
active force. This force guarantees phase separation and 
introduces surface-tension effects jl6|. 

This model has been applied with considerable success 
in measuring contact angles [TjJ and the effect of wall wet- 
tabilities, topography and micro-structure on drag reduc- 
tion of fluid flow through micro-channe ls 1181 . As an ex- 
tension of the Shan-Chen model, Benzipjj first derived 
an analytical expression for the contact angle and the 
surface energy between any two of the liquid, solid and 
vapor phases. 

In the Shan-Chen multiphase model, the non-ideal in- 
teraction is obtained by using an attractive short-range 
force 

F(x) = -Gty(x, t) ^ wM* + c*At, t)c u (9) 

i 

Adhesive forces between the fluid and solid phases are 
modeled by introducing an extra force 

F Q ds(x) = -G Qds V(x, t) ^ Wl s{yi + Ci Ai, t)c u (10) 

i 

where 
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Here s(x, t) = 0, 1 for the fluid and the solid phase, re- 
spectively, G, the interaction strength, is used to control 
the two-phase liquid interaction, and the adhesion pa- 
rameter G a d s is used to control the wettability behavior 
of the liquid at solid surfaces, ^(x^) is a local 'effective 
mass' [II 0, US 11, which is defined as[2l| 



v={r- ^)c 2 s At, 



(8) 



V*(x,t) = Voe- po/p , 



(12) 
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Using these definitions, the fluid momentum is changed 
at each time step according to 

pa =^c 4 /, + r(F + F Qds ), (13) 

i 

where u is the new fluid velocity. The equation of state 
in the Shan- Chen model is[2~i| 

P = p RT+^Mp)} 2 , (14) 

III. NUMERICAL SIMULATION OF CONTACT 
ANGLE 

The LB simulations were carried out in a 2D do- 
main. The grid mesh used is 50 x 200. In the simula- 
tion, the general non-slip bounce-back scheme was em- 
ployed for the solid-fluid interfaces, and periodic bound- 
ary conditions were applied at both inlet end and outlet 
end along the horizontal direction. A droplet, whose di- 
ameter is 30, is set on the middle between two ends. After 
30000 time steps (in our simulations the units are lattice 
units as Ref. [2l|, and the same as below), the result tends 
to stabilize. FIG. Q] shows different static contact angles 
for ^(x,i) = 4e- 20 °/'>[2l|, in Eqs. © and Values of 
parameters G and G a ds and the simulated contact angles 
(in degree) for each case are listed in TABLE fl] 




FIG. 1: Static contact angle. 



TABLE I: Adhesion parameters G a ds and the contact 
angles (in degree) for droplets. 



Case 


Adhesion 
parameters 


Contact angle 
for G = -120 


Contact angle 
for G = -130 


a 


-100 


147.2 


147.8 


b 


-130 


127.6 


132.4 


c 


-160 


109.8 


117.6 


d 


-190 


92.6 


102.6 





-220 


75.3 


87.6 


f 


-250 


58.6 


73.1 


g 


-280 


40.7 


59.5 


h 


-310 


18.0 


45.3 



As shown in FIG. [U G and G a ds determine the value 
of contact angle. G a ds represents the strength of interac- 
tion between fluid and solid surface and G represents the 
strength of interaction between fluid particles. A neg- 
ative G a ds indicates attractive interaction. When G is 
given, the greater the magnitude of |G a( j s |, the stronger 



the reaction, and thus resulting in smaller contact an- 
gle. So we can change parameter G and parameter G a d s 
to simulate contact angle for arbitrary surface condition, 
and then obtain different wall wettabilities. Form TA- 
BLE [l] we see that both parameters of G and G a ds have 
significant impact on the simulated contact angle. 

IV. NUMERICAL SIMULATION OF SLIP 
LENGTH 

In order to investigate the slip effect of boundary, we 
conducted numerical simulations of 2D Poiseuille flows. 
Typical density and velocity profiles of pressure-driven 
Poiseuille flows are displayed in FIG. [5] and FIG. [3] (given 
G = —120). The horizontal coordinate in both FIGs. 
represents the distance measured from one of the solid 
surface boundary to the other. The vertical coordinate 
in FIG. is the normalized velocity, where vo is the max 
velocity measured at the center in the channel for the case 
of no slip. The vertical coordinate in FIG. [3] is the nor- 
malized fluid density, where rhoo is the density of liquid 
for the case of no slip. The pressure gradient is specified 
as 0.005 in lattice unit for both cases. Different contact 
angles (as shown in both FIGs [2] and [3]) can be simulated 
by specifying different values of the adhesion parameter 
(G a( ; s ). All simulations were run until static equilibrium 
was nearly attained, and then a pressure gradient of 0.005 
was applied in the x-direction (flow direction). 




pi 

10 20 30 40 

y 



FIG. 2: Velocity profiles with different contact angles. 

As shown in FIG. [2j fluid velocity approaches zero as 
y — > (the lower boundary) or y — > 49 (the upper bound- 
ary) , which is consistent with the bounce-back boundary 
condition specified at the two boundaries. However, the 
fluid velocity increases dramatically in a very thin layer 
near the boundary. The layer is so thin that it is hard to 
see such details near the boundary, and velocity at the 
boundary looks like non-zero when plotted in a larger 
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given[l] by 




v{y) 



1 dP 

2/j, dx 



(h 2 -y 2 -2h6), 



(15) 



FIG. 3: Density profiles with different contact angles. 



scale, as shown by plots in FIG. H In a micro-scale, 
velocity at a boundary is zero, but in a macro-scale, the 
velocity appears non-zero (so called apparent slip). Plots 
shown in FIG. [5] clearly indicate that the slip velocity at 
the boundary increases as contact angle increases. In or- 
der to understand the physical mechanism of such kind of 
phenomenon, we also drew density profiles of fluid with 
different contact angles in FIG. [3] We should note that 
the density of the fluid with zero contact angle is con- 
stant (as shown by the horizontal line in FIG. How- 
ever, a sharp reduction of fluid density near the bound- 
ary is observed for a fluid with non-zero contact angle, 
which clearly indicates a layer of much less dense fluid 
(most probably gas) is induced between the dense liq- 
uid and solid surface. As discussed above, the parameter 
of G a ds controls the interaction between the fluid and 
solid surface. Increasing G a ds decreases the attraction 
(or increases the repulsion) between the liquid and solid 
surface, and thus attracts more gas to the surface. The 
less the dense of the fluid at the surface, the less viscous 
shear force, and the more significant slippage. 

From FIG. [5] and FIG. |3]we can see that the wetting 
properties of the wall have an important influence on the 
velocity profile, especially the slip velocity at the bound- 
ary. The more hydrophobic the wall is, the larger the 
slip velocity is on the boundary. There is a low-density 
layer between the bulk liquid and the wall, and the more 
hydrophobic the wall is, the lower density of fluid is, see 
FIG. [3J This result is similar to those obtained from 
other LB model simulation [22j and observed in MDS(23|. 
Compared with the macroscopic flow, the ratio of the 
low-density layer region to the inner region is larger in 
the microscopic flow, and this is the main difference be- 
tween micro-flow and macro-flow. Thus, the slip cannot 
be ignored in the micro-flow. 

Following Navier's hypothesis, the velocity in flow 
direction at position y between two parallel planes is 



where 2h is the distance between the planes, /i is the vis- 
cosity, dP/dx is pressure gradient, and S is slip length. 
Both slip velocity (v ) and slip length (5) can be esti- 
mated by matching the numerical LM simulation results 
of velocity profile with Eq. ([15)) . 

FIG. |3] shows that slip velocity increases linearly as 
pressure gradient increases. The larger the contact angle 
is, the larger the slip velocity is achieved, as shown by 
the solid diamonds (a contact angle of 127.6° obtained 
by using G = —120 and G a d s — —130) and squares (a 
contact angle of 58.6° obtained by using G = —120 and 
Gads = —250). However, the results shown in FIG. [5] 
clearly indicate that the estimated slip length is inde- 
pendent of pressure gradient. 
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FIG. 4: The slip velocity against pressure gradient. 
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FIG. 5: The slip length against pressure gradient. 



V. THE RELATIONSHIP BETWEEN SLIP 
LENGTH AND CONTACT ANGLE 

As discussed above, given an interaction strength (G), 
different contact angles and slip lengths can be simulated 
by using a series values of the adhesion parameter (G a ds) 
(see the values specified in TABLE H| . Numerical results 
of contact angles and slip lengths, represented by differ- 
ent symbols with different colors shown in FIG. [HI are cor- 
responding to different values of G, respectively, squares 
(G = -120), triangles (G = -125), triangles (G = -130) 
and triangles (G = —135). The numerical results shown 
in FIG. [HI cover a wide range of contact angles, ranging 
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from 18° to 150°. Though completely different interac- 
tion strength between fluid particles (G) and different 
interaction strength between fluid and solid {G a ds) were 
used to simulate both contact angle and slip length, our 
numerical results (as shown in FIG. [6]) clearly indicate 
that slip length is a function of contact angle. The rela- 
tionship between them can be expressed as Eq. (|16[) (as 
shown by the solid curve in FIG. which is obtained 
by fitting numerical results. 

5 L = 0.1441(e e/56 06 - 1) (16) 

where 9 is contact angle (in degree), and 5l is slip length 
in lattice unit. Our results are also consistent with those 
of Zhang and Kwok|22j (see solid circle in FIG.[B]for com- 
parison). As discussed above, the slip length is the result 




Contact angle 

FIG. 6: The slip length against contact angle. 

of interaction between the liquid and the wall surface, 
and it depends on the properties of the liquid (e.g. G) 
and the wall surface (e.g. G a ds)- So the simulation of 
slip length should take into account all parameters of the 
interaction objects. Though it is possible to construct 
a model of computing slip length in dependence of all 
interaction parameters, it is a very hard task [24|. In 
fact, the contact angle is a more comprehensive expres- 
sion of interactions between liquid and solid surface, and 
the weaker the solid-fluid interaction is, the larger the 
contact angle is. So we use contact angle as control pa- 
rameter to simulate the slip length. The result listed in 
FIG. [6] shows that this method is feasible, practical and 
much simpler. 

VI. THE UPPER LIMIT OF NANO- PARTICLES' 
DIAMETER FOR DRAG REDUCTION 

A practical application of Eq. fp~6|) is the determina- 
tion of the upper limit of nano-particles' size for the pur- 
pose of drag reduction. Water flooding has been widely 
used worldwide as a secondary recovery method, and has 
been proved an effective method. However, many techni- 
cal challenges have been encountered during field tests 
and applications. One of them is the extremely high 
injection pressure, especially for tight formations. Di. 
et.al.Q proposed an effective way of improving water in- 
jection and reducing drag of fluid flowing through rock's 



micro-channels with application of nano-particles adsorp- 
tion method, and have investigated the drag reduction 
mechanism experimentally and theoretically. When so- 
lution containing hydrophobic nano-particles (HNPs) of 
S1O2 is injected into a reservoir, HNPs are adsorbed to 
the wall of micro-channels to form a strong or super hy- 
drophobic layer, which can lead to a slip boundary con- 
dition and thus decrease drag. The slip model is shown 
in FIG. 

In FIG. [7aJ 2h represents the original pore throat di- 
ameter, and Iq is slip length brought about by substrate 
which is usually hydrophilic and has contact angle in 
the range of to 30 degree. When HNPs with effec- 
tive diameter of l p are adsorbed to the wall of the pore 
throat and formed a single-layer of nano-particles, the 
pore throat diameter decrease from 2h to 2(h — l p ), as 
shown in FIG. I7bl u s represents slip velocity in both 
FIG. |7al and FIG. EB 




(a) Original configuration (b) Configuration after HNPs 

adsorption treatment 



FIG. 7: The slip model. 

In general, l p does not equal the diameter of nano- 
particle because of nano-particle's intrinsic contact an- 
gle. After the treatment of HNPs adsorption, the wall's 
wettability is changed from hydrophilic to hydrophobic. 
The hydrophobicity of the solid surface increases the con- 
tact angle significantly and thus yields a much larger slip 
length 8. The necessary and sufficient condition for drag 
reduction can be represented as 

S>(!p + h), (17) 

For a 2D micro-channel, we assumed that nano- 
particles adsorbed on the wall are arranged closely in 
a single-layer, as shown in FIG. [8j where 6 p is the nano- 
particles' intrinsic contact angle. According to the prin- 
ciple of minimization of the Gibbs energy of a system [25j], 
it is easy to determine the location of liquid-gas surface, 
or the effective diameter of HNPs, l p , as shown in FIG.[8l 
lp = d p /2{l—cos6 p ), d p is the diameter of a nano-particle. 
Eq. (fT7|) becomes, 

5>^-(l-cos6 p ) + l , (18) 
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The upper limit of nano-particles' diameter for drag 
reduction can be obtained by rearranging Eq. ([T? 



2(5 -l ) 



^pmax 



(19) 



The slip length in Eq. ([T9]) , 8 , can be calculated with 




'///7///////7/////////////////7////////7//// 



FIG. 8: Schematic of nano-particles adsorbed on the 
wall. 

Eq. p^|) . but the contact angle of a rough surface (as 
shown in FIG. [8]), , is an unknown parameter. Accord- 
ing to the Cassie-Baxter equation [26|], the contact angle 
for the heterogeneous wetting regime is given as 



cosO = ficosd p - f 2 , 



(20) 



where f± and f% are the area fraction of the liquid-solid 
contact and liquid-gas contact, respectively, and can be 
determined by following equations. 



only when the diameter of nano particles is smaller than 
34 nm. 



We should note that both Eqs. (|T9J and J23J) are de- 
rived under the assumption that all nano-particles have 
same diameter and are well arranged as shown in FIG. [5] 
However, it is not the case for practical applications. In 
fact, the nano-particles have a statistical range of diam- 
eter, and they cannot be so well arranged and adsorbed 
on the wall. The wall surface is much more rough. So 
the above upper limit is a theoretical result and will be 
different from the practical situations. But according to 
Lauga, E., et al[28j, a more rough surface may lead to 
a greater contact angle. This conclusion has also been 
confirmed by our experiments, where the largest contact 
angle is 148.1°, as shown in FIG.O Therefore, Eq. 
yields a more conservative result. Though it is not an 
exact solution, the upper limit estimated by Eq. (1191) is 
still useful for determination the suitable size of nano- 
particles for a given reservoir. 




h = 7t - e p , (21) 

f 2 = 1 - sin6 p , (22) 

Then Eq. (|20|) can be expressed as 

cosO = (-7T — 9 p )cos9p + sinOp — 1, (23) 

In Eq. ([2"3[) . the contact angle 9 is independent of the 
diameter of nano-particle and only depends on nano- 
particle's intrinsic contact angle, when nano-particles ad- 
sorbed on the wall are arranged closely in a single-layer as 
shown in FIG. [H If the intrinsic contact angle 9 p equals 
to 120° ( is usually less than 120° [13] ), the apparent con- 
tact angle 9 equals 131.115° from Eq. ([23]). That is to 
say, when the nano-particles adsorbed on the wall are 
as shown in FIG. [51 the largest apparent contact angle 
may reach 131.115. Then a slip length of 1.35 lu (in our 
simulations the units are lattice units as Ref. [2lh can be 
got according to Eq. (TT5)) . In most cases, the substrate is 
hydrophilic and usually has contact angle in a range of 
to 30°and the corresponding slip length is smaller than 
0.102 lu. In our example, 1 lu= 20.408nm. When the in- 
trinsic contact angle of nano-particles is 120°, the upper 
limit of the diameter of nano-particles is about 34nm. 
Thus, the nano-particles adsorbing method is effective 



FIG. 9: The contact angle of a core sample after 
adsorbing nano-particles is 148.1 degree. 



VII. CONCLUSIONS 

(1) Both contact angle of a fluid on a solid surface and 
slip length of a fluid flowing over a solid surface can be 
easily simulated with the Shan-Chen LB model. 

(2) Both the non-ideal interaction between fluid parti- 
cles (controlled by a parameter of G) and the interaction 
between a fluid and a solid (controlled by another param- 
eter of G a ds) have significant impact on contact angle and 
slip length. However, they have little influence on the re- 
lationship between contact angle and slip length. 

(3) There is an upper limit on the diameter of nano- 
particles when applying nano-particle adsorption method 
to reduce drag. The upper limit can be estimated by 
applying the relationship between contact angle and slip 
length given the properties of nano-particles and geomet- 
rical parameters of a micro-channel. 
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